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Final  Technical  Report 

Pulse  Propagation  in  Random  Media  F  49620-95-1-0137 
95/02/01  -  98/01/31 

The  research  conducted  during  this  period  was  initially  done  in  collaboration  with 
George  Papanicolaou  (Stanford  University)  and  Benjamin  White  (Exxon  Research  & 
Engineering  Co.)  During  the  last  half  of  this  time  period,  however,  the  work  was  done  in 
collaboration  with  Benjamin  White  and  Leonard  Smka  (Exxon  Exploration  Co.).  In  broad 
terms,  the  work  focused  upon  direct  and  inverse  problems  involving  wave  propagation  in 
complex  heterogeneous  media,  modeled  as  random  media.  More  specifically,  Ae  work  has 
focused  upon  the  following  areas: 

(1)  Multiple  scattering  of  electromagnetic  waves  in  randomly  layered  dissipative  media.  The 
magnetotelluric  (MT)  problem  considered  has  lent  itself  to  modeling  within  the  “weak 
noise”  or  Gauss  Markov  regime.  In  addition  to  considering  propagation  per  se,  the  issue 
of  the  detection  of  anomalous  layers  has  been  addressed.  Details  are  presented  in  the 
enclosed  preprint,  “Random  Scattering  in  Magnetotellurics”,  by  Benjamin  White,  Werner 
Kohler  and  Leonard  Smka,  which  has  been  submitted  for  publication  to  Geophysics. 

(2)  Extension  of  the  ideas  mentioned  above  to  the  case  of  an  undulating  (locally  layered) 
medium.  The  theory  developed  for  layered  media  is  shown  to  be  robust. 

(3)  Incorporation  of  the  above  analysis  into  an  inversion  strategy. 

1.  Magnetotelluric  Probing 

We  consider  plane  electromagnetic  waves  normally  incident  upon  a  layered  Earth 
model  in  which  the  resistivity  is  assumed  to  vary  with  depth.  This  is  Ae  standard  MT 
model.  Very  low  frequencies  are  considered,  i.e.  10"^  <  /  <  10^ .  At  these  low  frequencies, 
skin  depth  is  large  and  the  incident  EM  energy  is  capable  of  penetrating  to  considerable 
depth.  Certain  substances,  such  as  groundwater  and  hydrocarbons,  have  anomalously  large 
resistivities  while  others,  such  as  hot,  pressurized  brine,  are  very  conducting.  One  of  the 
problems  we  considered  was  the  detection  of  such  anomalous  layers. 

The  main  contribution  of  this  work  is  our  accounting  for  the  effects  of  rapidly 
varying,  fine  scale  resistivity  variations.  The  problem  falls  within  the  theory  developed  by 
Khasminskii  ( Theory  of  Probability  and  its  Applications,  1 1,  1966,  p.  211-228).  Using 
this  theory,  the  surface  impedance  (measured  at  a  number  of  discrete  frequencies  across  the 
band  of  interest)  is  shown  to  be  asymptotically  described  as  a  Gaussian  random  vector.  The 
mean  value  of  this  vector  is  the  “efective  medium”  impedance  vector,  i.e  it  is  the 
impedance  vector  one  would  measure  were  the  randondy  fluctuating  resistivity  replaced  by 
a  smooth  effective  medium.  The  covariance  matrix  of  the  surface  impedance  fluctuations  is 
explicitly  computable.  The  theory  is  robust  and  does  not  depend  upon  any  detailed 
assumptions  about  the  resistivity  microstructure  fluctuation  statistics. 

The  Gaussian  character  of  the  limiting  impedance  statistics  makes  the  signal 
processing  aspects  of  the  problem  particularly  tractable.  For  the  problem  we  consider,  a 
thin  highly  resistive  layer  is  assumed  to  either  be  present  or  not  at  some  depth.  The  task  at 
hand  is  to  use  the  measured  surface  impedance  data  to  reach  a  decision.  The  Neyman- 
Pearson  lemma  provides  the  optimal  test.  As  was  mentioned,  the  Gaussian  character  of  the 
limiting  statistics  makes  this  test  relatively  straightforward  to  implement. 

Although  the  study  conducted  deds  exclusively  with  the  MT  problem,  the  ideas  are 
general  in  scope  and  would  apply  to  other  scenarios  in  which  the  underlying  scales  of  the 
problem  are  similar.  One  application  that  might  be  considered  is  the  problem  posed  by 


Richard  Albanese  several  years  ago,  in  which  the  detection  of  buried  chemical  pools  was 
the  issue. 


2.  Locally  Layered  Media: 

The  effort  here  has  been  to  extend  the  model  to  accomodate  deviations  from 
layering.  Undulations  are  introduced  and  the  goal  is  to  both  extend  the  relevant  theory  to 
this  case  and  also,  simultaneously,  to  demonstrate  the  robustness  of  the  layered  media 
results.  For  the  acoustic  problem  studied  at  length  in  SIAM  Review  (December  1991,  p. 
519-625),this  has  been  difficult;  a  number  of  false  starts  and  dead  ends  have  occurred. 
Moreover,  the  results  obtained  have  proven  too  compicated  to  be  useful. 

For  the  simpler  Gauss-Markov  limit  associated  with  the  above-mentioned  MT 
problem,  however,  the  situation  is  simpler.  Because  of  the  scaling,  a  straightforward 
perturbation  scheme  can  be  used,  as  shown  below,  to  account  for  the  randomness  and 
geometric  undulations  as  additive  effects.  In  fact ,  the  effects  of  the  randonmess  are 
ultimately  computed  using  the  plane-layered  model  since  the  effects  of  the  undulations  upon 
the  randomness  is  of  higher  order. 

Let  6  denote  a  small  parameter  as  in  the  preprint,  representing  the  ratio  of  the 
microscale  correlation  length  to  the  macroscale.  Assume  that  the  undulations  are  defined  as 

level  curves  ofz'  =  z  +  4dZ(x,y,z),  where  Z  is  a  known  function.  We  introduce  a  new 
coordinate  system  of  the  form: 

x'  =  x  +  4SX(x,y,z,4S),  y' =  y  +  ^fSY(x,y,z,4S),  z'  =  z  +  4SZ(x,y,z) 

where  X.Fmust  be  determined  so  that  x',y'are  everywhere  ±  z'.  Assuming  X,Y  to 

possess  expansions  in  powers  of  sfS ,  these  functions  can  be  determined.  (Generally,  the 
x',y'  coordinates  cannot  be  chosen  to  be  mutually  orthogonal,  however.)  ii  terms  of  these 
new  coordinates,  we  assume  that  the  conductivity  cr  is  a  rapidly  varying  function  of  z', 
i.e.  a  =  a{z'lS). 

We  can  then  derive  equations  for  the  electric  and  magnetic  fields  transverse  to  z'. 
The  derivations  are  tedious  but  straightforward  extensions  of  the  ideas  in  Marcuvitz’s 
Waveguide  Handbook.  We  obtain: 


B,E,.  =  -  S,.H,  +  S[x,3,Hy  -  -  Y„H,.)]  + 

Sa" (X.  -  - 3,.H,)-4SZ^E,.  - 4da-'zJd,H,  - 3,.H,.)  + 

+  Y,)h^.  +  (l  -  +  0(3) 


3,E,.  =  <y-'3,.[3,.H,.  -  +  4s(x,3,.Hy  -  Y,3,H^.  +  X^H,.  -  + 

■JS<t-'{y^-Z,)3,.{3,.H^ -3,.H,.)-^Z„E,-4s<t-%{3,.H^. -3,.H^.)- 
iaitiys{x,  +  +(i  -  ■Jsz,)h^.'\  +  0(S) 


E,  =  <y-'\3,H,  -  3,.H,.  +  4s{x,3,H,.  -  Y,3yH,  +  X^,  -  4^?,.)]+ 0(«) 


Analogous  equations  for  the  magnetic  iBeld  components  are  easily  obtained  by  duality. 

Stacking  the  (locally)  transverse  field  components  into  a  4  x  1  vector  X  leads  to  the 
following  equation  structure: 

y  1 


d,.X  = 


L,+46\^-^L,+L, 


X 


where  is  an  0(1)  deterministic  operator,  ^  A  ^  ^  appropriately  scaled  zero  mean 

random  operator  arising  from  the  random  resisitivity  fluctuations  and  L^isa.  deterministic 
operator  arising  because  of  the  undulations  in  the  layering.  To  leading  order,  the  effects  of 
the  undulations  occur  as  an  additive  operator.  Employing  a  perturbation  expansion  leads  to 
an  additive  contribution  due  to  the  undulations  at  the  level  of  the  fields  themselves.  The 


leading  order  term  corresponds  to  effective  medium  theory,  the  correction 


associated  with  A  ‘S  the  usual  Gauss  Markov  correction  associated  with  the  layered 


randomness  and  the  0\y6j  correction  associated  with  Lj  arises  due  to  the  undulations. 

This  simple  stracture  leads  to  explicitly  computable  results.  Computations  are 
currently  being  done  for  representative  geometries.  The  case  where  the  lateral  undulations 
are  themselves  random,  but  independent  of  and  slowly-varying  relative  to  the  rapid  random 
layering,  is  also  being  evaluated. When  complete,  the  results  will  be  submitted  for 
publication. 


3.  Inversion: 


The  results  discussed  in  (1),  and  presented  in  more  detail  in  the  enclosed  preprint, 
lend  themselves  naturally  to  the  problem  of  inversion.  In  this  case,  the  surface  impedance 
measurement  data  is  used  to  infer  the  large  scale  subsurface  resistivity  stmcture  (i.e.  the 
effective  medium  resistivity).  The  contribution  of  the  Gauss-Markov  theory  is  to  explicitly 
characterize  the  cross-correlations  of  the  measured  multi-frequency  surface  impedance 
noise. 

Consider  Maximum-Likelihood  estimation.  Because  of  the  multivariate  Gaussian 
character  of  the  measured  data,  the  task  reduces  to  minimizing  a  quadratic  functional  of  the 
form: 

(x,-X„)''C(x,-X,)  +  ln(detC.) 

where  Xj  is  the  measured  data  vector  while  and  C„  represent  the  mean  vector  and 
covariance  matrix  corresponding  to  the  particular  model  under  consideration.  In  particular, 
if  N  layers  (resting  upon  a  semi-infinite  basement)  are  being  used  in  the  model  to  attempt 

to  fit  the  data,  then  X^  and  C„  are  functions  of  3W+ 1  variables,  the  length,  mean 
conductivity  and  noise  strength  of  each  layer  together  with  the  basement  conductivity.  In 
Appendix  A  of  the  enclosed  preprint,  formulas  are  given  that  enable  one  to  evolve  the  mean 
impedance  as  well  as  the  covariance  matrix  across  a  layer.  The  explicit  nature  of  these 
formulas  enables  one  to  similarly  evolve  the  partial  derivatives  of  these  quantities  with 
respect  to  the  model  parameters.  In  this  way  one  obtains  an  explicit  evaluation  not  only  of 
the  objective  function  at  the  current  model  state  but  also  its  gradient  with  respect  to  the 
model  parameters. 

These  concepts  have  been  developed  and  are  currently  being  implemented.  When 
complete,  the  results  will  be  submitted  for  publication. 
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RANDOM  SCATTERING  IN  MAGNETOTELLURICS 
BENJAMIN  S.  WHITE,  WERNER  E.  KOHLER,  AND  LEONARD  J.  SRNKA 


Abstract.  Typical  well  logs  show  substantial  variations  of  formation  electrical  resistivity  over 
small  spatial  scales,  down  to  the  resolution  of  the  logging  tool.  Using  a  plane  stratified  Earth  model, 
we  examine  the  effects  of  this  fine  scale  microstructure  on  scattering  of  the  natmredly  occurring 
electromagnetic  waves  used  in  magnetotellurics  (MT).  We  show  theoretically  how  MT  data  may  be 
viewed  as  arising  largely  from  a  smoothed  “effective  medium”  version  of  the  resistivity  vs.  depth 
profile.  The  difference  between  the  data  produced  by  the  true  medium  and  that  produced  by  the 
effective  medium  is  due  to  random  scattering  noise,  which  is  fundamental  to  MT  since  it  arises 
from  the  very  small  spatial  scales  that  are  usually  ignored.  This  noise  is  substantial  and  hcis  unique 
statistical  properties,  which  we  characterize.  We  examine  the  implications  of  the  existence  of  this 
noise  for  the  detectability  of  a  thin  layer  of  Increased  resistivity  at  depth.  The  theory  is  shown  to 
agree  well  with  Monte  Carlo  simulations. 


1.  Introduction 

In  magnetotelluric  (MT)  surveys  [1],  [2],  [3],  surface  measurements  of  the  Earth’s  electrical 
impedance  over  a  broad  frequency  range  at  a  number  of  different  sites  are  analyzed  to  produce 
maps  of  electrical  resistivity  in  the  subsurface.  Naturally  occurring  ambient  electromagnetic  (EM) 
radiation  is  used  as  a  source.  The  method  is  widely  used  for  the  mapping  of  crustal-scale  structures 
both  on  land  and  at  the  seafloor,  and  can  also  be  useful  for  sedimentary  basin  exploration  in  areas 
where  seismic  reflection  data  quality  is  poor. 

The  principal  drawback  of  MT  is  its  limited  accuracy  and  spatial  resolution,  due  to  the  diffusive 
character  of  low  frequency  EM  propagation  in  the  Earth.  Mathematically  this  is  to  be  expected, 
since  the  processing  of  the  MT  surface  measurements  to  produce  a  subsurface  map  is  an  inverse 
problem  which  is  classically  ill-posed  in  the  sense  of  Hadamard  [4].  That  is,  many  different  solutions 
can  fit  the  observed  data  almost  as  well.  In  fact,  EM  prospecting  was  a  seminal  problem  in  inverse 
theory.  Early  research  on  EM  prospecting  was  done  by  Tikhonov  [5],  [6],  [7],  [8]  and  Cagniard 
[9];  the  theoretical  aspects  of  the  inverse  problem  were  later  studied  by  Weidelt  [10].  Tikhonov,  in 
particular  [5],  showed  that,  contrary  to  prevailing  wisdom,  ill-posed  problems  could  be  solved  both 
theoretically  and  in  practice  (for  a  history,  see  [2],  Chapter  1).  Today,  the  Tikhonov  regularization 
technique  has  achieved  great  success  in  solving  inverse  problems  in  many  diverse  scientific  fields 
[11]. 

The  basic  assumption  of  regularization  theory  is  that  the  solution  is  piecewise  smooth.  In 
choosing  among  candidate  solutions,  smoothness  is  valued  in  addition  to  the  candidate’s  ability  to 
explain  the  data,  and  it  is  this  additional  criterion  that  allows  a  “best”  solution  to  be  found.  In 
the  words  of  Parker  [12]  (p.  294),  “we  seek  the  simplest,  least  exciting  kind  of  solution”. 
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However,  the  smoothness  assumption  is  clearly  false  for  MT,  as  is  evident  from  even  a  casual 
inspection  of  well  logs.  As  an  example,  Figure  1  shows  formation  resistivity  as  a  function  of  depth 
for  an  onshore  Louisiana  oil  exploration  well  in  the  Gulf  of  Mexico  basin.  The  lithology  consists 
primarily  of  intermixed  layers  of  sandstone  and  shale.  Note  that  the  geologic  section  is  quite 
conductive,  which  is  typical  for  this  basin.  The  data  were  recorded  using  an  Atlas  Wireline  DIFL 
induction  logging  tool  (RILD  response),  with  a  nominal  two  foot  vertical  resolution.  The  figure 
illustrates  the  weU  known  fact  that  formation  resistivity  is  rapidly  varying  on  very  small  spatial 
scales  down  to  the  spatial  resolution  of  the  tool.  Furthermore,  the  amplitude  of  the  resistivity 
fiuctuations  is  by  no  means  small.  In  general,  resistivities  varying  rapidly  over  one  or  more  orders 
of  magnitude  are  not  uncommon  in  well  log  data.  In  this  paper  we  examine  the  effects  of  this  fine 
scale  structure  on  EM  scattering,  and  the  subsequent  implications  for  the  fundamental  limits  of 
vertical  resolution  in  MT  data  processing. 


O.^  0.6  O-S  i  -1-2  1-4.  1.6  1.S  a 
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Figure  1.  Onshore  Louisiana  Well  Log  Data:  Resistivity  (Ohm-m)  vs.  Depth  (ft). 

The  record  extends  from  3,500  ft  to  19,900  ft  in  depth. 

It  is  natural  to  describe  the  small  scale  resistivity  fluctuations  as  a  random  process  called  the 
microstructure.  An  MT  survey  will  not  recover  resistivity  with  all  the  resolution  of  the  microstruc¬ 
ture,  but  only  some  kind  of  smoothed,  spatially  averaged,  version.  This  idea  will  be  made  more 
precise  below,  where  fturthermore  we  show  that  it  is  conductivity,  i.e.  the  reciprocal  of  resistivity, 
that  is  averaged.  But  to  explore  the  spatial  scales  involved,  consider  a  sliding  window  with  a  length 
of  500  ft  centered  at  different  depths  in  Figure  1.  As  the  center  of  the  window  moves  along  the 
abscissa  in  Figure  1,  the  average  resistivity  over  that  window  changes  smoothly.  Besides  the  aver¬ 
age,  other  estimated  statistics  will  also  vary  with  window  position.  As  a  typical  example.  Figure  2 
shows  the  normalized  resistivity  autocorrelation  function,  estimated  from  a  500  ft  window  centered 
at  a  depth  of  4, 450  ft.  From  this  figure,  the  length  scale  of  the  microstructure  can  be  estimated 
as  about  1  =  13  ft  (nominally  4  m).  This  is  the  microscale,  i.e.,  the  correlation  length  of  the 
microstructure;  two  points  separated  by  much  more  that  4  m  will  have  resistivity  fluctuations  that 
are  not  significantly  correlated,  while  significant  correlations  will  exist  for  distances  much  smaller 
than  this.  Our  studies  of  other  portions  of  this  log,  and  of  other  resistivity  logs  with  different 
lithologies,  have  produced  qualitatively  similar  results. 
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0  to  20  30  40  so  60 

Figure  2.  Typical  Onshore  Louisiana  Well  Log  Autocorrelation  Function:  Nor¬ 
malized  Resistivity  Autocorrelation  vs.  Lag  (ft),  measured  over  a  500  ft  window 
centered  at  a  depth  of  4,450  ft. 


This  value  of  I  is  consistent  with  studies  of  acoustic  well  logs,  which,  while  recording  a  different 
material  property,  may  be  expected  to  vary  on  the  same  spatial  scale.  In  reflection  seismology, 
scattering  from  thin  beds  is  explained  by  the  now  classical  theory  of  O’Doherty  and  Anstey  [13]. 
Their  Figure  11  shows  an  autocorrelation  function  of  reflectivities  that  is  at  least  qualitatively 
consistent  with  an  exponential  autocorrelation  function  of  acoustic  impedance,  and  a  correlation 
length  I  of  a  few  meters.  Godfrey  et.  al.  [14]  and  Banik  et.  al.  [15]  use  decreasing  exponentials  to 
get  good  data  fits  to  the  autocorrelation  function,  and  give  explanations  for  the  observed  expo¬ 
nential  behavior  in  terms  of  Markovian  models  for  the  sedimentation  process.  Banik  et.  al.  cite 
a  typical  exponential  decay  length,  corresponding  to  our  f,  of  4.5  m.  White  et.  al.  [16]  also  foimd 
autocorrelation  functions  that  were  well  fit  by  exponentials,  with  f  varying  between  1  m  and  10  m, 
and  a  most  typical  value  of  about  3  m.  Other  studies  of  the  statistics  of  acoustic  well  logs  that  are 
generally  consistent  with  these  results  are  reported  by  Walden  and  Hosken  [17],  Velzeboer  [18],  and 
Sato  [19].  Some  of  the  larger  values  of  I  obtained  by  these  authors  have  been  up  to  20  m.  A  review 
of  the  mathematics  of  acoustic  scattering  in  stratified  media  with  fine  scale  structure  is  given  by 
Asch  et.  al.  [20]. 

In  what  follows,  we  consider  horizontally  plane  stratified  media,  so  that  resistivity  varies  as  a 
function  of  depth  only.  Besides  the  microscale  /,  we  consider  a  macroscale  L,  on  which  significant 
variations  of  resistivity  might  be  recovered  by  MT.  The  resistivity  is  then  modeled  as  a  function  of 
these  two  disparate  spatial  scales  L»l.  The  ratio  5  =  l/L  is  a  small  parameter  which  is  used  in 
the  asymptotic  analysis  that  follows. 

In  Section  2  we  show  under  wide  hypotheses  that,  despite  the  microstructure,  there  is  a  piecewise 
smooth  resistivity  vs.  depth  curve  that  approximately  explains  the  data.  More  specifically,  we  use 
a  theorem  of  Khasminskii  [21]  to  derive  “effective  medium  theory”  —  i.e.  we  construct,  by  suitable 
averaging,  a  smooth,  effective  medium  which  will  produce  surface  impedance  measurements  that 
match  those  produced  by  the  true  resistivity  profile,  but  with  a  small  error  of  order  0{V6).  In 
this  context,  the  ill-posedness  of  the  inverse  problem  can  be  seen  as  a  consequence  of  effective 
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medium  theory.  Since  many  different  microstructures  can  be  smoothed  to  produce  the  same  effective 
medium,  it  is  not  possible  to  distinguish  between  these  different  microstructures,  as  a  complete 
inversion  algorithm  would  have  to  do.  The  inverse  problem  becomes  well-posed  if  one  solves  for 
the  effective  medium  only.  However  in  this  case  the  data  cannot  match  the  solution  of  the  forward 
problem,  but  will  be  in  error  by  an  amount  of  order  0{y/6). 

This  order  0{VS)  error  in  the  data  is  a  kind  of  “geologic  noise”  caused  by  multiple  scattering 
from  the  fine  scale  microstructure.  This  noise  is  fimdamental  to  MT  and  cannot  be  reduced  by 
instrumentation,  since  it  is  produced  solely  by  a  discrepancy  between  the  true  resistivity  vs.  depth 
profile  and  the  smooth  curve  that  is  assumed  in  the  mathematical  inversion.  Moreover,  this  kind 
of  scattering  noise  can  be  substantial.  In  Section  4  (See  Figure  3)  we  compare  an  estimate  of  the 
magnitude  of  geologic  noise  to  that  of  all  other  sources  of  error,  based  on  an  analysis  of  Parker  [12], 
and  find  that,  except  for  very  low  frequencies,  the  magnitudes  are  comparable. 

In  Section  2  we  also  apply  another  theorem  of  Khasminskii’s  [21]  to  characterize  the  unique 
statistical  properties  of  geologic  scattering  noise.  Further  formulas  for  the  effective  calculation  of 
these  statistics  are  developed  in  Section  3.  These  statistics  do  not  depend  on  the  details  of  the 
microstructrue,  but  only  on  some  averaged  quantities.  Furthermore,  because  Khasminskii’s  theorem 
is  a  kind  of  Central  Limit  Theorem  for  differential  equations,  the  siuface  impedance'scattering  noise 
has  a  Gaussian  distribution  even  when  the  resistivity  fluctuations  do  not.  Thus  there  is  a  good 
deal  of  model  independence,  and  it  is  not  necessary  to  undertake  a  detailed  statistical  analysis  of 
the  Earth’s  resistivity  microstructure  in  order  to  apply  the  theory.  In  particular,  for  the  case  of  a 
homogeneous  half  space  containing  random  fluctuations  {i.e.  a  random  half  space),  we  can  derive 
explicit  algebraic  formulas  that  predict  the  variances  and  covariances  of  the  real  and  imaginary 
parts  of  an  array  of  smrface  impedances  measured  at  different  firequencies.  These  formulas  depend 
only  on  a  single  microstructure  parameter  which  characterizes  the  noise  strength  and  merely  scales 
all  the  variances  and  covariances.  Thus,  for  a  homogeneous  random  half  space,  the  statistics  are 
universal. 

In  Section  4  we  study  the  implications  of  the  theory  for  the  detectability  of  a  thin  layer  of  material 
in  the  subsurface,  when  the  material’s  resistivity  varies  markedly  from  that  of  the  background,  as 
in  prospecting  for  hydrocarbons.  Ultimately,  detectability  is  also  a  question  of  spatial  resolution 
for  MT,  since  very  large  amounts  of  the  material  will  almost  certainly  be  detected,  while  very 
small  amounts  will  almost  certainly  be  missed.  Clearly,  detectability  also  depends  on  the  depth  of 
investigation. 

The  question  then  is,  “At  a  given  depth,  how  strong  a  target  is  necessary  for  detection,  given 
that  the  EM  waves  used  as  a  probe  are  scattered  by  the  microstructure  ?”  In  Section  4  we  develop 
the  theory  for  a  model  problem  of  this  type.  We  show  that  there  is  an  optimal  method  of  data 
processing,  based  on  the  Neyman-Pearson  lemma  [22],  [23].  Attempts  at  detection  will  make  two 
types  of  errors:  false  positives  (“detecting”  a  nonexistent  target)  and  false  negatives  (failing  to 
detect  a  target).  These  two  types  of  errors  cannot  be  eliminated,  but  can  be  traded  off  against 
each  other  by  setting  a  threshold  on  how  compelhng  the  evidence  needs  to  be  for  detection  to  be 
declared.  We  derive  expressions  for  this  tradeoff,  i.e.  expressions  relating  the  false  negative  rate  to 
the  false  positive  rate.  As  is  standard  practice  in  detection  theory  [23],  this  relation  between  error 
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rates  can  be  expressed  as  a  ROC  (Receiver  Operating  Characteristic)  curve,  which  is  a  plot  of  the 
probability  of  successful  detection  vs.  the  false  positive  rate.  Because  the  method  is  theoretically 
optimal,  its  ROC  curve  cannot  be  improved,  and  so  provides  a  fundamental  limit  to  detectability. 

In  Section  5,  we  compare  the  theory  to  Monte  Carlo  simulations.  We  obtain  excellent  agreement 
of  the  simulations  and  the  theoretical  noise  statistics  over  a  frequency  range  spanning  five  orders 
of  magnitude.  For  the  detectability  theory  we  include  both  geologic  scattering  noise  and  a  white 
noise  component,  to  model  other  sources  of  error.  The  detectability  theory  also  compares  well  with 
simulations. 

Although  we  have  not  pursued  the  idea  in  this  paper,  knowledge  of  the  noise  statistics  should  be 
of  direct  value  in  the  inversion  of  MT  data,  by  suggesting  the  theoretically  optimal  misfit  function 
to  be  used  in  optimization.  A  simple  example,  corresponding  to  current  practice,  is  the  use  of 
variance- weighted  least  squares  [12],  [24].  This  misfit  function  is  appropriate  if  the  data  errors 
at  each  frequency  are  uncorrelated,  but  have  diflFerent  variances.  However,  for  geologic  scattering 
noise,  the  errors  at  difierent  frequencies  are  highly  correlated,  as  we  show  here.  So  the  appropriate 
misfit  function  should  be  a  quadratic  form  in  the  errors,  weighted  by  all  the  elements  of  the  inverse 
of  the  covariance  matrix.  Whether  substantial  improvement  in  inversions  can  be  obtained  using 
this  misfit  function  is  a  subject  for  future  research. 

2.  Spatial  Scales  and  a  Stochastic  Limit 

We  consider  a  plane  stratified  Earth  occupying  the  half  space  z  <  0,  with  a  basement  layer 
consisting  of  a  homogeneous  half  space  for  z  <  —L  and  with  air  in  z  >  0.  Let  a  be  conductivity, 
H  the  magnetic  permeability,  e  the  dielectric  permittivity,  E  the  electric  field  and  H  the  magnetic 
field,  fj,  will  be  assumed  to  be  constant  while  a  and  e  will  be  assumed  to  vary  with  z,  but  not  with 
frequency.  For  angular  frequency  w  and  time  dependence  Maxwell’s  equations  become 

V  X  E  =  iw/uH 

V  X  H  =  (or  -  twejE.  (2.1) 

For  a  normally  incident  plane  wave  (with  s-directed  electric  field  and  y-directed  magnetic  field), 
we  obtain  the  usual  equations  for  the  scalar  components  E  and  H 

—iunH 

(cr  —  ioje)E.  (2.2) 

To  nondimensionalize,  we  choose  typical  constants  /i,  i,  u,  a,  H  and  define  a  typical  length  scale 
L,  impedance  ^  and  reference  electric  field  strength  E  by 

L  =  l/y/uiiia 

i  =  l/aL 

E  =  in.  (2.3) 

Note  that  y/2L  is  the  skin  depth.  Let  the  impedance  be 

^  =  E/H. 


(2.4) 
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Magnetotelluric  surveys  produce  measiuements  of  ^  at  z  =  0  over  a  range  of  frequencies  a;. 

Set 

z'  =  z/L,  h'  —  plIJx,  e' —  e/li  u]'  =  u/(jj 

H'  =  H/H,  E'  =  EIE,  ^  =  ^/|.  (2.5) 

Then  equations  (2.2)  and  (2.4)  hold  again  for  the  primed  variables. 

In  what  follows,  we  use  the  nondimensional  variables  but,  for  notational  convenience,  drop  the 
primes.  Also,  for  z  <  0,  the  parameter  regimes  to  be  considered  will  permit  us  to  employ  the  usual 
quasi-static  approximation  [2];  that  is,  displacement  currents  are  neglected  so  that  effectively  e  may 
be  set  equal  to  0.  Then,  substitution  of  (2.4)  into  (2.2)  yields 


— ^  -  iwju, 


z  <  0.  (2.6) 

Because  of  the  continuity  of  E  and  H,  ^  is  also  continuous  across  layer  boundaries.  The  condition 


that  there  are  no  upgoing  waves  in  z  <  —L  yields  that 

E  =  ^bH  at  z  — L, 

where 


(2.7) 


ks 


m 


(2.8) 


= 

kjB  = 

Thus,  an  initial  condition  for  (2.6)  is 

^  at  z  =  -L.  (2.9) 

We  taJce  permeability  to  be  constant  and  consider  variations  in  the  conductivity.  In  particular, 
we  consider  deterministic  macroscale  variations  of  the  conductivity  on  a  spatial  scale  of  order  i, 
and  random  microscale  variations  on  a  much  smaller  spatial  scale  I  «  L.  Let 

S  =  l/L  «  1,  (2.10) 

so  that  in  dimensional  units,  a  =  a  {zll^zji^.  Then  in  nondimensional  units,  z'  =  zjL.  We 
obtain,  after  again  dropping  primes  for  notational  convenience,  that 

a  =  (j{z/5,z) .  (2-11) 

The  random  function  (t(zi,Z2)  now  drives  equation  (2.6)  stochastically,  with  z\  =  z/6,  Z2  =  z. 
We  assume  that,  for  each  fixed  Z2,  a  (zi,Z2)  is  a  stochastic  process  in  zi  which  satisfies  a  mixing 
condition.  This  means  that  cr{zi  4-  Azi,Z2)  becomes  statistically  independent  of  o'(zi,Z2)  as  the 
microscale  separation  distance  Azi  oo.  In  other  words,  we  assume  that  the  microstructure 
has  no  long-range  correlations.  Other  technical  conditions  are  given  by  Khasminskii  [21],  whose 
theorems  we  will  now  apply. 

Let^i,^25---  ,  be  iV  impedances  corresponding,  respectively,  to  the  JV  frequencies  a;i,a;2>  ••  •  ,<^N- 
We  will  compute  their  joint  statistics,  approximately,  for  small  J.  Let 

(’i  =  9fe}.  j  = 


(2.12) 
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We  define  the  vector  X  e  by 

x  =  (2.13) 

Then,  equation  (2.6)  for  the  complex  impedance  ^  implies  that  X  satisfies  a  vector  stochastic 
differential  equation  of  the  form 

£x  =  F{z/6,  X),  X(-L)  =  Xo-  (2.14) 

Here 


Xo  —  (^B,l)  •  •  •  •  •  •  >^B,N>  (2.15) 

where  are  the  real  and  imaginary  parts  of  ^b-,  given  by  (2.8),  for  frequency  Wj,  and 

F(s,X)  =  <7(s,X2n+i)H(X)+V 
H(X)  =  (((?f -((If, 

V  =  (0 . 0,  -uiM, -unit,  if  ■  (2.16) 


Equation  (2.14)  is  of  the  type  considered  by  Khasminskii  [21].  To  apply  his  Theorem  (1.1),  we 
compute 

F(X)=  hm  i  r£?{F(s,X)}ds,  (2.17) 

z-^oo  z  jQ 

where  E{-}  denotes  expected  or  mean  value.  The  deterministic  function  F(X)  is  the  average,  both 
in  probability  and  space  (over  the  microstructure  only),  of  F(s,X).  Then,  by  Theorem  (1.1)  X 
may  be  approximated,  to  leading  order  as  5  j.  0,  by  the  nonrandom  vector  X,  which  satisfies 

=  F(X) 

X(-L)  =  Xo-  (2.18) 


Equation  (2.18)  is  equivalent  to  effective  medium  theory.  To  see  this,  note  from  (2.16)  that  F  is 
obtained  from  F  by  replacing  the  random  process  cr{s,  z)  by  its  local  average  over  the  microstructure, 
d{z),  where 


cr{z)  =  lim  -  f  E{a{s' ,  z)}ds' . 
S-+00  s  Jq 


(2.19) 


Therefore,  each  is  approximated,  to  leading  order,  by  the  deterministic  impedance  which 
satisfies 


-iujjfi,  -L<z<0 

^j{-L)  =  ^B,j  (2.20) 


Note  that  d,  and  hence  vary  only  on  the  macroscale  since  the  local  averaging  (2.19)  has  erased 
all  microscale  variation. 

To  derive  corrections  to  effective  medium  theory,  we  next  apply  Khasminskii’s  Theorem  (3.1). 
Let 
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Then  as  5  |  0,  the  statistics  of  Y  are  well-approximated  by  the  statistics  of  Y*^, 
Gaussian  random  process  that  satisfies  the  linear  integral  equation 

where  Y°  is  a 

Y»(z)  =  W»(z)  +  Ap  (X(f))  Y<i(^)dC. 

(2.22) 

Here  W°(z)  is  a  Gaussian  process  with  independent  increments,  zero  expectation  and  a  correlation 
matrix 

f;{(W«(z))  (W»(z))’'}  =  A  (X(0)  dC, 

(2.23) 

where 

A(X)=  lim  -  r  rE{[P(si,X)-F(X)]  [P  (s2,X)  -  F(X)]^}dsids2, 

Z-+OO  Z  Jq  Jq 

(2.24) 

Prom  (2.16),  (2.24)  we  obtain  that 

A  (X)  =  7"  (X2N+1) ,  H  (X)  H  (X)"^ , 

(2.25) 

where 

j^{z)  =  lim  -  /  /  E{{a  (si,  z)  -  d{z))  {a  {s2,  z)  -  d{z))}dsids2. 

S-^00  S  Jq  Jq 

(2.26) 

A  more  convenient  representation  for  the  correction  term  is  obtained  by  expressing  W°(z)  as  an 
Ito  stochastic  integral  [25] 

w°(z)  =  y‘%(C)H(X(C))d/?(C), 

(2.27) 

where  I3{z)  is  a  standard  Brownian  motion  (the  Wiener  process).  Using  the  “white  noise” -like 
property  of  the  Brownian  diflferential 


E{dm)dm)}  =  <^(Ci  -  C2)rfCi,  (2.28) 

(where  J(Ci  —  C2)  is  the  standard  (J-function).  We  may  verify  that  W°  defined  by  (2.27)  has  the 
required  covariance  (2.23). 

Using  (2.27),  equation  (2.22)  can  be  differentiated  to  give  a  set  of  Ito  stochastic  differential 
equations,  i.e.  “white  noise”  equations.  We  identify  the  corrections  as 

(2.29) 

so  that,  from  (2.21)  written  in  component  form,  we  obtain 

^3  (^)  =  !?■  (^)  +  li  {z)  •  (2.30) 

The  equations  for  the  are 

=  -‘^d{z)lj{z)^j{z)dz+'i{z)^j{z)dl3{z) 
ij{-L)  =  0.  (2.31) 

Thus,  the  fy  are  linear  “white  noise”  perturbations  of  (2.20).  In  fact,  (2.31)  can  be  derived  heuris- 
tically  by  substituting  (2.30)  into  (2.6),  assuming  that  satisfies  the  effective  medium  equations 
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(2.20).  Then,  one  obtains,  to  order  0{^/S) 

(2.32) 

Equation  (2.31)  is  obtained  from  (2.32)  by  replacing  the  first  (t{z15,z)  by  its  average  and  by 
identifying  [a{zl5,z)  —  d’(z)]/v/5  as  an  appropriate  white  noise  with  covariance  proportional  to 
7^(z),  which  is  the  “efiective  noise  strength”  of  the  microstructure  at  macroscopic  position  z. 

To  summarize,  the  conductivity  is  first  averaged  locally,  via  (2.19),  and  the  result  is  used  to 
compute  the  impedance  in  the  effective  medium  of  conductivity  ct(z),  via  equation  (2.20).  The 
leading  order  correction,  obtained  via  (2.30),  is  of  order  0{y/5),  and  is  given  by  the  which 
satisfy  the  Ito  stochastic  differential  equations  (2.31).  We  shall  next  use  (2.31)  to  compute  the 
joint  statistics  of  the 

3.  Random  Scattering  Statistics:  Propagation  Equations  and  a  Homogeneous 

Random  Half-Space 

Prom  equation  (2.31),  it  follows  that  the  are  zero  mean,  jointly  Gaussian  random  variables, 
and  so  it  remains  to  compute  their  covariances  at  z  =  0  to  determine  the  statistics  of  the  impedance 
fluctuations  at  the  surface  due  to  multiple  scattering  from  the  microstructiure.  Let 

cfi^  =  = 

C7f  =  =  j,l  =  l,2...,iV.  (3.1) 

Let  *  denote  complex  conjugate.  We  wifi  determine  the  real  matrices  in  (3.1)  from  the  two  complex 
matrices 

Cji=E{ijti}-  (3.2) 

Note  the  symmetry  properties;  C  is  complex  symmetric  and  C  is  Hermitian. 

=  C,  =  C.  (3.3) 

Then,  from  the  identity  -1-  iij  we  obtain 

=  isR{C  +  C},  =  isRlC  -  C} 

=  ^SS{C-C},  =  i9{C -f- C}.  (3.4) 

To  derive  an  equation  for  C,  we  use  Ito’s  multiplication  table  [25],  Table  1. 


dp 

dz 

d/3 

dz 

0 

dz 

0 

0 

Table  1.  Ito’s  Multiplication  Table 


Then 


d  (ijii)  =  ildij  +  ijdii  +  dijdii- 


(3.5) 
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In  the  stochastic  calculus,  the  extra  term  d^jd^i  is  not  zero  since  {d(3)^  =  dz  according  to  Table  1. 
Substituting  from  (2.30)  into  (3.5)  we  obtain 

d  (I, -6)  =  -2^  {ij  +  ii)  ijiidz  +  7  (lig  +  61?)  dP  +  'T^^i^dz.  (3.6) 

Now,  taking  expected  values  in  (3.6),  and  using  E{dP}  =  0,  we  obtain 

^Cjiiz)  =  -2a(6  +  ei)C'j7+7"§l? 

Cjli-L)  =  0.  (3.7) 

The  initial  condition  in  (3.7)  follows  from  (2.31),  i.e.  the  absence  of  random  fluctuations  at  z  =  —L. 
Similarly,  from  the  stochastic  difierential 

d  {ijit)  =  itdij + ijdit  +  {di^m)  (3.8) 

we  obtain 

Cji{-L)  =  0.  (3.9) 

As  an  important  special  case,  consider  a  homogeneous  random  half  space,  z  <  0.  Then  a  and 
7  are  constants  and  we  can  let  fr  ->  oo.  In  this  case,  the  initial  conditions  are  ignored  and  the 
z-independent  solutions  for  C  and  C  are  found  from  (3.7)  and  (3.9)  to  be 

„2\  /  \  /^2\ 


=(£) 


■'ji 


{±' 

V2S, 


,  0  C 


(3.10) 


To  determine  we  may  set  d^j/dz  =  0  in  (2.20)  to  get  an  algebraic  relation  for  we  obtain 


thorn  _  “i7r/4 

“V  ^ 

Substitution  of  (3.11)  into  (3.10)  yields 

/  ^2^3/2 


(3.11) 


^hom  _ 


23/2^5/2 
23/2^5/2  j 

Substitution  of  (3.12)  into  (3.4)  yields 


UjOJl 


\  (-1  -  i) 


rfhom  _ 


-  v^)] 

^  (jjj  +  (jjl 


)■ 


(3.12) 


QRR{hom)  _  aji-^/u}jU}i,  _  Q.ji[u}j  +  U)1  +  y/UJjUli) 


where 


C,RIlhom) 
jl  = 


aji 


^IR(hom) 

Cji  =  -Oj/Wi, 


/  ] 

/  OJjUl 

123/2^/2  J 

V  (wj  +  wO(v^ + \M)  / 

(3.13) 


(3.14) 


Traditionally,  magnetotelluric  data  has  been  displayed  in  terms  of  the  apparent  resistivity  [1] 

UJfJl  ’ 


pa 


(3.15) 


where  ^  is  evaluated  at  the  surface  z  =  0. 
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Substitution  of  (2.29)  into  (3.15)  yields 

Pa  «  Pa  +  pa, 


where 


(3.16) 


(3.17) 


Pa  = 


UIJL 


(3.18) 


From  (3.18)  we  can  compute  the  covariance  of  pa  at  two  different  frequencies,  u)j  and  coi,  as 

=  (;;^)  miSCj, + CjiiCj,}.  (3.19) 


To  evaluate  (3.17),  (3.19)  for  a  homogeneous  random  half  space,  put  (3.11)  into  (3.17)  to  obtain 

Pa"^  =  4.  (3.20) 

CT 

That  is,  pa°”*  is  independent  of  frequency  and  equals  the  resistivity  of  the  effective  medium.  How¬ 
ever,  the  random  fluctuations  in  apparent  resistivity  are  not  independent  of  frequency.  Substitution 
of  (3.12)  into  (3.17)  yields 


ipf  -/lom 

MPa,j  Pa,l  )-7 


,  IpUjUJl  ( 

V  2^7  V( 


1 


+ 


+  ■y/Sj')  ^ 


xy^ + v^)  + ^i) 

In  particular,  for  Wj  =  wj  =  w  we  obtain  the  variance  of 


(3.21) 


(3.22) 


The  variance  increases  as  the  square  root  of  frequency,  and  so  is  zero  at  a;  =  0.  This  can  be 
understood  as  resulting  from  an  averaging  over  lengths  of  the  order  of  a  skin  depth 


=  P- 

V  wpcr 


(3.23) 


As  frequency  decreases,  skin  depth  increases  and  the  wave  averages  over  more  microstructure. 
Consequently,  the  variance  decreases;  effective  medium  theory  becomes  more  accurate.  Inserting 

(3.23)  into  (3.22)  yields 


(3.24) 


The  variance  is  thus  inversely  proportional  to  skin  depth  and  proportional  to  noise  strength  7^. 
The  5'“'*  dependence  then  follows  from  dimensional  analysis. 

The  homogeneous  random  half  space  enables  us  to  study  the  effects  of  random  conductivity 
fluctuations  in  a  setting  free  of  any  macroscale  artifacts.  How  important  are  these  fluctuations  as 
a  source  of  error? 

Parker  [12],  in  Table  5.02A  of  his  text,  presents  an  MT  data  set  that  includes  a  tabulation  of 
both  apparent  resistivity  and  a  standard  deviation  of  uncertainty  over  a  broad  range  of  periods. 
We’ll  use  this  data  to  provide  a  benchmark  for  an  order-of-magnitude  comparison. 
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Prom  (3.20)  and  (3.24)  we  form  the  relative  error  (the  standard  deviation  divided  by  the  mean) 


^  T^hom 

ra 


(3.25) 


In  order  to  evaluate  this  expression,  we  consider  a  model  consisting  of  layers  of  thickness  1.  The 
conductivities  of  these  layers  are  assumed  to  be  independent,  identically  distributed  random  vari¬ 
ables  uniformly  distributed  on  the  interval  [<y\,(T^  (where  1,  ai  and  a2  represent  actual,  dimensional 
values).  For  this  model,  Ra  reduces  to 


Ra  — 


(3.26) 


where  Zs  likewise  represents  the  actual  dimensional  skin  depth  for  the  effective  medium  half  space. 

Figure  3  compares  Ra  with  corresponding  apparent  resistivity  quotients  formed  from  Parker’s 
data  and  plotted  as  a  function  of  logarithmic  frequency.  We  have  used  the  values  o‘2  =  0.1  S/m, 
(Ti  =  0.01  S/m  and  f  =  3  m. 


Figure  3.  Comparison  of  Uncertainty  Levels:  Apparent  Resistivity  Uncertainty 
vs.  Logarithm  of  Frequency  (Hz).  The  solid  curve  represents  uncertainty  due  to 
resistivity  microstructure.  The  circles  are  uncertainty  levels  reported  by  Parker  [12]. 

The  sources  of  uncertainty  that  we  compare  in  Figure  3  are  certainly  different.  However,  this 
graphical  comparison  does  illustrate  the  fact  that  a  reasonable  level  of  conductivity  fluctuations 
can  give  rise  to  uncertainties  that  are  roughly  comparable  in  size  over  a  significant  frequency  band 
to  those  encountered  in  practice. 

4.  Fundamental  Limits  to  Detectability 

In  the  interpretation  of  magnetotelluric  surveys  it  is  often  assumed  that  the  Earth  consists  of  a 
discrete  number  of  plane  layers,  within  each  of  which  the  conductivity  cr  is  a  constant.  This  type 
of  model  ignores  small  spatial  scale  random  variations,  which  are  always  present.  In  this  section 
we  investigate  the  effect  of  scattering  noise  due  to  the  microstructure  on  the  ability  of  such  a 
survey  to  detect  a  particular  layer  of  interest.  Since  scattering  noise  is  created  by  the  same  physical 
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mechanism  as  the  “signal”  from  the  macroscopic  layers,  this  type  of  noise  represents  a  fimdamental 
limit  to  detectability  that  cannot  be  reduced  by  instrumentation. 

For  a  more  general  plane-layered  model  that  includes  the  microstructure,  we  consider  cr  to  be  a 
piecewise-homogeneous  rcindom  function.  We  consider  m  macroscopic  layers,  with  layer  boundaries 
—L  =  Zm  <  Zm-i  <  Zm-2  <  . . .  <  zq  =  0.  For  Zj  <  z  <  %-i, (7  =  cTj{z/S),  where  aj  is  a  stationary 
random  process  with  mean  aj  and  noise  strength  7?.  Formulas  for  rapid  calculation  of  the  effective 
medium  impedances  and  of  the  covariances  for  models  of  this  type  are  derived  in  Appendix  A. 

For  simplicity,  we  consider  in  detail  the  case  m  =  3,  illustrated  by  Figure  4.  The  top  layer,  zi  < 
z  <  zo  has  an  average  conductivity  di  and  an  effective  scattering  strength  7|.  These  same  values 
are  assumed  to  be  taken  by  the  third  layer,  which  lies  directly  above  the  basement.  Sandwiched 
between  these  two  layers,  in  ^2  <  ■^  <  ^1,  is  a  layer  with  different  values,  ^2,72. 
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Figure  4.  Random  Three-Layer  Model:  Resistivity  (Ohm-m)  vs.  Depth  (m).  A 
1000  Ohm-m  layer  (100  m  thick)  lies  between  two  layers  having  randomly-fluctuating 
microstructure.  The  semi-infinite  basement  has  a  1000  Ohm-m  resistivity. 


We  investigate  here  the  statistical  reliability  with  which  the  middle  layer  can  be  detected  with 
an  MT  survey.  The  existence  of  this  layer,  j.e.  the  model  shown  in  Figure  4,  is  referred  to  as 
hypothesis  Hi.  This  hypothesis  is  to  be  contrasted  with  the  one  in  which  the  middle  layer  is  not 
present,  i.e.  there  is  just  one  layer,  occupying  <  z  <  zq,  with  parameters  di  and  71-  This  latter 
hypothesis  is  referred  to  as  Hq,  the  null  hypothesis. 

This  example  is  a  model  for  MT  exploration.  Since  hot,  pressurized  brine  is  highly  conducting, 
d2  >  di  in  prospecting  for  geothermal  energy.  Since  hydrocarbons  are  highly  resistive,  ^2  <  di  in 
prospecting  for  oil  or  gas.  In  either  case,  we  are  to  reach  a  decision  Hi  or  Hq,  as  to  whether  or  not 
the  target  has  been  detected  by  MT. 

Within  this  framework,  there  is  an  optimal  method  of  reaching  a  decision  by  processing  the 
survey  data.  This  “most  powerful  test”  is  given  by  the  Lemma  of  Neyman  and  Pearson  [22].  Let 
the  data  vector  be 

d  =  i^n)  •  (4.1) 

As  before,  are,  respectively,  the  real  and  imaginary  parts  of  the  impedance  at  frequencies  Uj 
measured  at  the  smrface  z  =  0. 
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For  /i  =  0  or  1,  let  denote  the  expectation  under  the  assumption  of  hypothesis  Hh-  Then,  for 
each  case,  the  mean  of  the  data  is 

,|N,hill,h)l2,h)---  )lN,h)  ■  (4.2) 

That  is,  the  2N  vector  of  real  and  imaginary  parts  of  the  impedances  for  the  effective  medium 
corresponding  to  the  appropriate  hypothesis. 

Similarly,  putting  the  covariance  matrix  for  the  data  depends  upon  the  hypothesis 

{Hi  or  Ho).  For  each  hypothesis  we  have 


(/^RR  /^RI  \ 

cf  J  • 


(4.3) 


That  is,  the  2N  x  2N  covariance  matrix  is  partitioned  into  four  N  x  N  blocks,  which  can  be 
computed  for  the  appropriate  model  using  the  general  theory. 

Since  the  data  are  Gaussian,  the  probability  density  of  the  data  vector  in  21V’-dimensional  space 
is,  for  /i  =  0  or  1, 

exp[-4(d  -  dh)'^Cj7^(d  -  dh)] 


Phid)  = 


(4.4) 


■v/det27rCh 

The  “most  powerful  test”  involves  the  choice  of  a  threshold  9.  Given  this  threshold,  we  are  to 
decide  Hi  whenever 

Pi(d) 


Po(d) 


>9. 


(4.5) 


For  the  Gaussian  of  (4.4),  this  test  is  equivalent  to  thresholding  the  quadratic  form  Q(d),  i.e. 

Q(d)  =  (d-do)'^Co'(d-do)-(d-di)'^Ci-'(d-di)>0.  (4.6) 

Here,  the  alternative  threshold  9  for  Q  is  related  to  9  by  the  relation 


0  =  2  log  0  +  log  det  Cl  —  log  det  Cq.  (4.7) 

Now,  (3(d)  has  a  probability  density  qo{Q)  if  hypothesis  Hq  is  true  and  a  probability  density 
qiiQ)  if  hypothesis  Hi  is  true.  The  false  positive  rate  Ppp  is  the  probability  of  deciding  that  a 
target  layer  has  been  detected  when,  in  fact,  it  is  not  present.  From  (4.6) 

Ppp  =  Pfp{9)  =  f  qo{Q)dQ.  (4.8) 

Je 

The  false  negative  rate,  Pfn,  is  the  probability  that  the  decision  made  will  miss  a  target  that  is 
present;  so,  from  (4.6) 

Pfn  =  Pfn{&)  ~  f  qiiQ)dQ.  (4.9) 

J  — oo 

These  two  types  of  errors  can  be  traded  off  against  each  other  by  changing  the  threshold  9.  In 
general,  however,  one  error  rate  cannot  be  reduced  without  increasing  the  other  error  rate.  The 
trade-off  relation  between  false  positive  emd  false  negative  error  rates  is  a  fundamental  limitation 
on  the  reliability  of  detection.  It  is  standard  practice  in  detection  theory  to  illustrate  the  relation 
between  the  error  rates  by  plotting  a  ROC  (Receiver  Operating  Characteristic)  curve,  in  which  the 
probability  of  detection,  1  —  Pfn,  is  plotted  against  the  false  positive  rate,  Pfp  [23].  Since  the 
Neyman-Peaxson  Lemma  guarantees  that  our  data  processing  is  optimal,  the  ROC  curve  cannot  be 
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improved.  These  ideas  are  illustrated  in  the  next  section,  where  the  theory  is  compared  to  Monte 
Carlo  simulations. 

The  above  formulas  simplify  considerably  if  we  make  the  approximation 


Cl  «  Co  =  C.  (4.10) 

This  is  the  case,  for  instance,  when  the  middle  (target)  layer  is  small,  and/or  the  contrast  in 
properties  is  not  great,  i.e.  when  detection  is  difficult.  If  (4.10)  is  valid,  then  the  quadratic  terms 
in  (4.6)  cancel  and  Q  becomes  linear,  i.e. 

Q  «  2(di-do)^C-M-d7c-^di+dJC-^do 

=  (di  -  do)^C-H2d  -  di  -  do).  (4.11) 


Under  hypotheses  Hq  and  Hi,  respectively,  linearized  Q  has  mean  values 


Qo  =  £;o{Q}  =  -(di-do)^C-^(di-do) 

Qi  =  Ei{Q}  =  (di  -  do)^C-ndi  -  do).  (4.12) 

Under  either  hypothesis,  Q  has  standard  deviation  S,  where 

S  =  2y^(di  -  do)^C-i(di  -  do)  =  2y/^i.  (4.13) 


Using  (4.11)  and  (4.12),  the  false  positive  and  false  negative  error  rates  for  the  optimal  test  (4.6) 
can  be  represented  using  the  cumulative  normal  probability  distribution  function 


/: 


exp(— T^/2)dr. 


(4.14) 


Then 


Ppp  =  Pfp{0)  =  \—F  ^  — 

PpN  =  Pfn{0)  =  F  .  (4.15) 

Approximation  (4.10)  also  simplifies  the  analysis  of  stacked  data.  Suppose  that  data  vectors 
of  the  form  (4.1)  are  collected  firom  r  different  locations,  widely  separated  on  the  surface.  If  it 
is  assumed  that  the  macrostructure  is  identical  at  each  location,  but  that  the  microstructures  are 
statistically  independent  (with  identical  probabilities),  then  the  data  vectors  may  be  considered  to 
be  r  independent,  identically  distributed  samples.  The  optimal  test  then  utilizes  only  the  average 
of  the  r  vectors;  this  average  is  used  in  place  of  d  in  equation  (4.11).  The  resulting  Q  is  then 
compared  to  a  threshold.  Stacking  reduces  the  errors;  the  stacked  error  rates  become  ‘ 


ppp^’^{e) 

pp^^{e) 


i-f  (.(.?--9o)y^) 


(4.16) 
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5.  Monte  Carlo  Simulations 

We  use  the  three-layer  model  shown  in  Figure  4  as  a  prototype  to  numerically  study  and  evaluate 
the  theory.  A  semi-infinite  basement  having  a  constant  1000  Ohm-m  resistivity  is  assumed  to 
occupy  the  region  z  <  —6000  m.  A  thin,  highly  resistive  layer  of  interest  (100  m  thick  with  1000 
Ohm-m  resistivity)  is  positioned  at  a  depth  of  2000  m,  i.e.  —2100  m  <  0  <  —2000  m.  The 
intervening  regions,  —6000  m  <  ^  <  -2100  m  and  —2000  m  <  z  <  0  m  consist  of  random  layers 
having  independent,  identically  distributed  conductivities.  Since  MT  is  quite  insensitive  to  resistive 
targets  in  conductive  backgrounds,  this  is  a  good  test.  Recall  that  the  theory  possesses  a  robustness 
with  respect  to  microstructure  details;  it  predicts  that  the  statistical  behavior  of  the  impedance  and 
apparent  resistivity  depends  only  upon  the  mean  conductivity  a{z)  and  the  fluctuation  parameter 
S'y'^  and  not  upon  the  details  of  the  statistical  model  generating  these  parameters.  To  provide 
some  illustration  of  this  robustness,  we  consider  two  statistical  models.  In  the  first  model,  the  layer 
conductivities  are  uniformly  distributed  on  the  interval  0.01  S/m<z<  0.10  S/m  and  the  layers  are 
each  assumed  to  be  a  constant  3  m  in  thickness.  In  the  second  model,  the  layer  thicknesses  are  also 
randomized  with  a  uniform  distribution  on  the  interval  1  m  <  /  <  5  m.  The  layer  conductivities 
are  again  assumed  to  be  independent  and  identically,  uniformly  distributed.  In  this  case,  since  the 
mean  layer  thickness  is  again  3  m,  one  can  show  that  the  conductivity  interval  generating  the  same 
values  of  a  and  remains  the  same. 

Two  cases,  where  the  thin  resistive  layer  is  both  present  and  absent,  are  considered.  (When  the 
thin  layer  is  absent,  the  random  layering  is  simply  taken  to  occupy  the  entire  interval  —6000  m  <  z  < 
0  m.)  5000  realizations  were  generated  for  each  case  {i.e.  thin  layer  present  and  absent)  and  for  each 
random  layering  model.  Fifty  one  fi-equencies,  equally  spaced  on  a  logarithmic  scale  extending  firom 
—2  to  3,  were  considered.  The  results  are  displayed  in  Figure  5.  We  plot  the  apparent  resistivity 
mean  and  the  mean  ±  one  standard  deviation  over  the  frequency  range  10“^  <  /  <  10^.  Simulation 
results  for  the  two  random  models  are  compared  with  the  theoretical  predictions.  Agreement  is 
good  over  the  entire  five  decade  firequency  band.  As  one  would  expect,  the  resistivity  fluctuations 
decrease  as  frequency  decreases.  Moreover,  as  frequency  decreases,  the  entire  6  kilometers  of  random 
layering  (with  the  thin  resistive  target  layer  present  or  not)  becomes  increasingly  transparent  and 
the  surface  apparent  resistivity  tends  toward  the  1000  Ohm-m  value  of  the  semi-infinite  basement. 

Figures  6  ofier  additional  evidence  of  good  agreement  between  the  theory  and  simulations  for 
the  real  and  imaginary  parts  of  the  impedance  and  for  the  apparent  resistivity.  The  first  two 
plots  compare  predicted  (Gaussian)  and  computed  probability  density  functions  for  the  real  cind 
imaginary  parts  of  the  impedance,  respectively,  at  1  Hz.  Note  that  for  both  impedance  components 
the  variance  is  small,  on  the  order  of  10“^,  so  that  the  peak  density  value  is  correspondingly  on 
the  order  of  thousands.  The  imaginary  part  is  negative  and  roughly  equal  in  magnitude  to  its 
real  counterpart,  as  one  would  expect  (i.e.  the  impedance  phase  is  roughly  —45°).  The  third  of 
Figures  6  presents  an  analogous  comparison  of  theory  and  simulations  for  the  apparent  resistivity 
at  1  Hz.  (For  simplicity,  the  simulations  used  in  these  and  all  subsequent  plots  are  those  involving 
the  fixed  3  m  layer  thickness.)  ' 
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Figure  5.  Apparent  Resistivity  (Ohm-m)  vs.  Logarithm  of  Frequency  (Hz):  a  com¬ 
parison  of  theory  and  simulations.  The  mean  and  mean  ±  one  standard  deviation 
are  plotted  for  (a)  the  thin  resistive  layer  present  and  (b)  the  layer  absent.  The  solid 
lines  represent  the  theory.  The  circles  and  stars  represent  the  two  random  layering 
models  used  in  the  simulations. 


Figure  6.  Impedance  Statistics  at  IHz:  Comparison  of  Theory  and  Simulations. 
Probability  density  functions  predicted  by  the  theory  (solid  curves)  and  constructed 
from  the  simulations  (circles)  are  shown  for  (a)  the  real  part  of  the  impedance,  (b) 
the  imaginary  part  of  the  impedance  and  (c)  the  apparent  resisitivity. 


Since  the  number  of  frequencies  iV  =  51  in  our  calculations,  covariance  matricies  Ch,  h  =  0, 1, 
as  given  by  (4.3),  are  102  x  102  real  symmetric  cirrays.  In  Figures  7  we  display  surface  plots  of  the 
covariance  matrix  subblocks  Cf®-,  and  as  functions  of  logarithmic  frequency.  Although 
the  particular  case  shown  corresponds  to  that  computed  from  the  simulations  with  the  thin,  highly 
resistive  layer  present,  the  other  possibilities  (theory  vs.  simulations,  layer  present  vs.  layer  absent) 
all  produce  indistinguishable  results  on  the  scale  shown.  The  behavior  shown  is  what  one  would 
expect.  As  frequency  increases,  random  fluctuations  in  both  the  real  and  imaginary  parts  of  the 
impedance  increase  in  amplitude.  Therefore,  as  one  moves  out  along  the  main  diagonal  in  the  first 
and  third  plots,  the  corresponding  variances  likewise  increase;  variances  of  the  imaginary  part  are 
larger  than  their  real  counterparts.  As  the  second  of  Figures  7  shows,  the  real-imaginary  cross¬ 
correlations  are  likewise  largest  in  amplitude  at  the  highest  frequencies.  Note  that  the  behavior 
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of  the  fluctuations  is  coupled  so  that  this  fluctuation  product  is  negative.  (This  negative  value 
agrees  with  analogous  cross-correlation  results  obtained  for  the  homogeneous  random  half  space; 
c.f.  (3.13),  (3.14). 


Figure  7.  Covariance  Submatrices:  (a)  Real-Real  (b)  Real-Imaginary  C^, 
and  (c)  Imaginary-Imaginary  Covariance  Submatrices  are  plotted  against  Log¬ 
arithm  of  Frequency. 


Figures  7  suggest  a  large  spread  in  the  amplitudes  of  the  covariance  matrix  elements  but  mask 
the  low  frequency  fine  structure.  Such  a  large  range  of  amplitudes  is  not  surprising,  given  the  results 
of  Section  3  for  the  homogeneous  random  half  space.  Note  in  particular  the  frequency  dependence 
in  (3.13)  and  (3.14),  together  with  the  fact  that  the  frequency  range  considered  spans  five  decades. 
The  matrix  elements  in  Figures  7  have  amplitudes  ranging  roughly  from  6  x  10"“*  to  6  x  10“^^. 

Figures  8  study  the  real-real  covariance  subblock  as  a  representative  case;  results  for  the 
other  subblocks  are  qualitatively  similar.  The  first  two  plots  display  the  low  frequency  structure 
for  theory  and  simulations,  respectively.  The  third  and  fourth  plots,  on  the  other  hand,  study 
the  impact  of  the  thin  resistive  layer  upon  this  covariance  subblock;  is  plotted  for 

both  theory  and  simulations.  These  plots  further  confirm  the  good  agreement  between  theory  and 
simulations.  They  also  lend  credence  to  approximation  (4.10) 

Low  frequency  smrface  impedance  information  is  very  important  in  detecting  the  presence  or  ab¬ 
sence  of  a  thin,  highly  resistive  layer  such  as  that  shown  in  Figure  4.  Higher  frequencies  correspond 
to  both  shorter  skin  depths  and  increased  impedance  fluctuation  levels  due  to  the  random  conduc¬ 
tivity  microstructure.  At  low  frequencies,  sources  of  noise  other  than  conductivity  microstructure 
assume  increasing  importance;  this  is  certainly  suggested  by  the  Parker  data  [12]  in  Figure  3. 

We  will  model  other  sources  of  noise  with  an  additional  white  noise  component.  That  is,  the 
measured  impedance  will  be  taken  to  be  . 

^j' +  ^ +  i{^j +^j) ,  j  =  (5.1) 

where  is  a  family  of  independent,  identically-distributed,  zero  mean  Gaussian  random 

variables  having  common  (dimensional)  variance  As  a  result  of  this  modification,  the  mean 
or  eff'ective  impedances  remain  unchanged  but  the  covariance  matrix  undergoes  the  change 
JC  5C  4-  uH. 

The  detection  problem  can  be  summarized  by  examining  Figure  9.  Assuming  =  10“®,  we 
have  plotted  the  mean  apparent  resistivity  and  the  mean  ±  one  standard  deviation  for  each  of  the 
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Figure  8.  Real-Real  Covariance  Submatrbc  vs.  Logarithmic  Frequency.  The  first 
two  figures  show  the  fine  scale  structure  of  as  computed  by  (a)  theory  emd 
(b)  simulations.  The  latter  two  figures  show  the  effect  of  the  thin  resistive  layer 
upon  the  covariance  submatrix,  t.e.  Cf^  —  as  computed  by  theory  in  (c)  and 
simulations  in  (d). 


two  cases;  the  solid-line  curves  represent  the  case  where  the  thin  resistive  layer  is  present  {Hi  true) 
while  the  dashed  curves  correspond  to  the  layer  being  absent  {Hq  true).  The  circles  and  stars,  on 
the  other  hand,  represent  two  sets  of  measurements,  in  this  case  two  simulation  realizations,  both 
with  the  layer  present.  Given  just  one  of  these  sets  of  measurements,  we  must  decide  whether  or 
not  the  layer  is  present. 

.  We  attack  this  detection  problem  using  the  ideas  of  Section  4.  In  order  to  determine  the  false 
positive  and  false  negative  rates  associated  with  this  problem,  i.e.  Ppp  and  PpN  as  given  by  (4.8), 
(4.9),  we  must  determine  the  probability  density  functions  and  qi  present  in  the  integrands. 
We  have  numerically  determined  these  density  fimctions  for  the  quadratic  functional  Q(d)  in  (4.6) 
and  found  that  the  approximation  Cq  w  Ci,  leading  to  (4.11),  is  a  very  good  one  for  our  model 
problem.  As  noted  in  Section  4,  this  linearizing  approximation  reduces  Q(d)  to  a  Gaussian  random 
variable  with  known  mean  and  variance  and  makes  the  theory  very  easy  to  implement.  Ultimately, 
we  eliminate  the  threshold  parameter  to  obtain  Ppp  vs.  Pfn.  This  relation  is  plotted  as  a  ROC 
curve,  i.e.  the  detection  probability,  1  —  PpN,  is  plotted, against  the  false  positive  rate,  Ppp. 
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Figure  9.  The  Detection  Problem:  Apparent  Resistivity  (Ohm-m)  vs.  Logarithm 
of  Frequency  (Hz).  The  solid  and  dashed  cmwes  are  predictions  of  the  theory  for  the 
resistive  layer  being  present  and  absent,  respectively.  The  mean  and  mean  ±  one 
standard  deviation  are  shown  for  each  case.  The  circles  and  stars  are  two  different 
realizations,  i.  e.  two  different  sets  of  random  microstructure,  both  for  the  case  where 
the  layer  is  present. 

To  illustrate  these  ideas,  consider  a  case  where  an  added  white  noise  level  =  10~®  exists 
and  the  highly  resistive  layer,  with  its  upper  surface  remaining  at  a  depth  of  2000  m,  has  its 
thickness  reduced  from  an  initial  100  m  to  a  final  10  m.  We  expect  the  layer  detection  to  become 
increasingly  difficult  as  its  thickness  is  reduced.  Figures  10  illustrate  these  ideas  at  the  level  of  the 
probability  density  functions  qh{Q)-  The  two  plots  correspond  to  layer  thicknesses  of  100  m  and 
25  m,  respectively.  In  the  first  case,  the  density  functions  qi  and  q^,  (corresponding  to  the  100  m 
layer  being  present  and  absent),  are  quite  separated.  Given  a  single  realization,  i.e.  a  single  set 
of  surface  impedance  measurements,  one  could  therefore,  with  relatively  high  probability,  correctly 
decide  upon  the  presence  or  absence  of  the  resistive  layer.  (If  the  density  functions  had  completely 
disjoint  supports,  a  single  realization  would  with  certainty  answer  the  detection  question.)  As  the 
layer  thickness  decreases,  the  density  functions  tend  to  peak  and  migrate  toward  each  other;  for 
the  25  m  thick  layer  shown  in  the  second  plot,  the  two  density  functions  have  considerable  overlap. 
In  this  case,  given  a  single  realization,  identifying  the  correct  hypothesis  is  correspondingly  more 
difficult.  In  the  extreme  case  where  both  density  functions  completely  overlap,  the  realization  {i.e. 
the  set  of  impedance  measurements)  would  provide  no  useful  detection  information. 

Symmetries  are  evident  in  Figures  10.  The  density  functions  go  and  gi  are  Gaussian  with  equal 
variances  and  means  symmetrically  located  relative  to  the  origin.  This  structure  foUows  from 
(4.11)-(4.13). 

Recall  that  Figure  9  includes  apparent  resistivity  plots  for  two  realizations  (i.e.  different  random 
microstructures),  both  with  the  100  m  resistive  layer  present.  One  can  ask,  in  the  context  of  the 
optimal  detection  theory  presented,  “How  detectable  is  the  resistive  layer  in  each  of  these  two 
particular  cases?”.  The  quadratic  form  Q(d)  defined  by  (4.6)  has  values  of  17.7  and  —2.1  for  the 
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Figure  10.  Probability  Density  Functions  for  the  Detection  Problem:  Density 
functions  go  and  Qi  for  a  (a)  100  m  thick  resistive  layer  and  (b)  25  m  layer.  In 
each  case,  the  density  function  qo{Q),  corresponding  to  the  layer  absent,  appears  on 
the  left  while  the  density  function  qi{Q),  corresponding  to  the  layer  present,  appears 
on  the  right. 


“circles”  and  “stars”  data,  respectively.  Locating  these  abscissa  values  on  the  first  of  Figures  10 
shows  that  the  17.7  value  is  well  within  the  probability  mass  of  qi,  essentially  disjoint  firom  the  the 
probability  mass  of  qo,  while  the  —2.1  value  is  in  the  overlap  region  of  the  go  and  gi  density  tails.  In 
the  former  case,  detection  of  the  resistive  layer  would  be  relatively  straightforward.  Noting  (4.9), 
6  =  17.7  corresponds  to  PpN  =  0.78;  any  false  negative  threshold  less  than  this  value  would  result 
in  the  correct  detection  of  the  layer.  On  the  other  hand,  6  =  —2.1  corresponds  to  PpN  =  0.01  and 
any  higher  threshold  would  cause  the  layer  to  be  missed. 

In  Figures  11,  the  performance  of  the  optimal  detector  is  illustrated  by  its  ROC  curve:  that  is, 
the  probability  of  detection,  1  —  Pfnj  is  plotted  against  the  false  positive  rate  Ppp-  When  the 
two  probability  densities  go  and  gi  are  well  separated,  good  detection  is  possible.  In  this  case  the 
probability  of  detection  is  close  to  one  unless  the  false  positive  rate  is  kept  very  small.  Therefore 
the  ROC  curve,  which  always  passes  through  the  origin,  rises  steeply  to  unity  in  this  case.  At  the 
other  extreme,  when  the  two  probabihty  densities  go  and  gi  substantially  overlap,  detection  will 
be  poor.  In  this  case  the  ROC  curve  Ues  close  to  the  diagonal  line  where  the  detection  probability 
equals  the  false  positive  rate;  hence  the  detector  performs  no  better  than  the  flip  of  a  biased  coin. 
Figures  11  clearly  show  a  progression  from  one  extreme  to  the  other  as  the  resistive  layer  thickness 
is  reduced  firom  100  m  to  10  m. 

One  can  similarly  study  the  detection  question  as  other  parameters  are  varied,  e.g.  one  where 
resistive  layer  thickness  is  held  fixed  but  the  level  of  additive  noise  is  varied.  Figure  12  shows  what 
happens,  for  the  case  of  the  50  m  layer,  as  the  level  of  additive  noise  is  increased  firom  =  10“^^ 
to  =  10“®.  At  the  lowest  10“^^  noise  level,  the  ROC  curve  predicts  almost  certain  detection, 
unless  the  false  positive  rate  is  kept  very  low.  At  the  highest  10“^  noise  level,  the  detection  curve 
is  basically  the  straight  line  of  equal  probabilities.  In  this  case  the  added  noise  has  corrupted  the 
measurements  to  the  point  that  they  no  longer  provide  useful  detection  information. 
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Figure  11.  ROC  Curves:  Probability  of  Detection  vs.  False  Positive  Rate  for  a 
resistive  layer  thickness  of  (a)  100  m,  (b)  50  m,  (c)  25  m,  (d)  10  m,  constructed 
using  the  theory.  The  circles  represent  a  corresponding  curve  created  using  the 
simulations. 
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Figure  12.  The  Effect  of  Added  Noise  upon  Detection:  ROC  curves  are  plotted  for 
a  resistive  layer  thickness  of  50  m  and  added  noise  levels  of  =  10“^^,  10“®,  10“^, 
and  10“®.  The  curves  decrease  monotonically  as  the  noise  level  is  increased. 


A  final  point  worth  emphasizing  is  the  important  role  that  the  correct  covariance  matrices  play 
in  the  detection  process.  Figure  13  illustrates  this  for  the  case  of  a  50  m  thick  resistive  layer 
and  an  additive  noise  level  of  =  10“^^.  The  solid  curve  is  the  ROC  curve  computed  using  the 
correct  covariance  matrix.  The  dcished  curve,  on  the  other  hand,  represents  the  performance  of  a 
detector  which  assumes  that  the  covariance  matrix  is  the  identity  matrix.  (Any  multiple  of  the 
identity  matrix  would  correspondingly  change  the  threshold  level  but  not  change  the  resulting  ROC 
curve.)  The  improvement  in  detection  capability  achieved  by  using  the  correct  covariance  matrix 
is  pronotmced. 
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Figure  13.  The  Importance  of  Using  Correct  Covariance  Information:  the  case  of 
detecting  a  50  m  thick  resistive  layer  in  the  presence  of  —  10“^^  added  noise  is 
considered.  The  solid  curve  is  the  ROC  curve  using  correct  covariance  information 
while  the  dashed  curve  represents  the  case  where  these  covariance  weights  have  been 
replaced  by  the  identity  matrix. 


6.  Conclusions 

The  Monte  Carlo  simulations  verify  the  accuracy  of  both  the  eflFective  medium  theory  and  the 
theory  of  random  scattering  statistics  over  a  frequency  range  spanning  five  orders  of  magnitude. 
Although  the  scattering  noise  is  significant  compared  with  estimates  of  other  noise  sources,  its  effects 
can  be  substantially  mitigated,  in  the  models  considered  here,  using  the  theory  of  this  paper.  If 
other  sources  of  noise  can  be  sufficiently  reduced,  good  detection  is  possible,  in  these  models,  for 
relatively  thin  layers  of  anomalous  resistivity  at  depth.  However,  good  detection  levels  are  achieved 
only  through  use  of  optimal  detection  algorithms,  which  incorporate  a  theoretical  understanding 
of  the  noise  statistics. 
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Appendix  A.  Piecewise-Homogeneous  Random  Media 

In  this  appendix  we  consider  piecewise-homogeneous  random  media  of  the  type  introduced  in 
Section  4.  These  types  of  models  generalize  the  usual  macroscopic  layer  models  to  include  the 
effects  of  random  scattering  from  the  microstructure. 

Consider  m  macroscopic  layers,  with  layer  boundaries  —L  =  Zm<  Zm-i  <  Zm-2  <  •  •  •  <  zq  =  0. 
For  Zr  <  z  <  Zr-i,(T  =  (Tr{z/5),  where  <Tr  is  a  stationary  random  process  with  mean  ar  and  noise 
strength  7^.  Consider  the  interface  at  z  =  for  some  r.  We  assume  that  ^{zr),C{zr)  and  C(^) 
have  been  determined.  (Recall  that  C  and  C  are  all  continuous  across  the  interface.)  We  will  next 
derive  formulas  for  |'(z),  C(z)  and  C(z)  for  Zr  <  z  <  z^-i-  Using  these  formulas  recursively,  one  can 
compute  the  relevant  quantities  at  z  =  0  rapidly,  by  stepping  up  from  the  basement  layer-by-layer. 
We  write  a  =  aril  =  Jr  and  z  =  Zr-  Let 


Vj^expfa  [  ^j(s)ds')  , 

(A.1) 

\  Jzo  / 

where  the  j-subscript  refers  to  frequency  ujj.  Then 

(A.2) 

and 

^Vi+kjVi  =  0, 

(A.3) 

where 

kj  =  y/ujua 

(A.4) 

Since  from  (A.l),  (A.2) 

II  II 

(A.5) 

we  obtain  from  (A.3),  (A.5)  that 


=  Aj  exp[iA:j(z  —  z)]  -I-  Bj  exp[—ikj{z  —  z)]. 


(A.6) 


where 


=  \( 

2  V‘  '  ikj 


(A.7) 


Solving  (A.2)  for  and  substituting  from  (A.6)  yields 

r  ^  /  Aj  exp[ikj{z  -  z)]  -  Bj  exp[-%(z  -  z)]\ 

^  a  \Ajexp[ikj{z  -  z)]+ Bjexp[-ikj{z  -  z)]J  '  '  ’ 

Prom  (A.8)  ^j  can  be  determined  at  the  next  layer  boundary,  z  =  Zr-i.  Formula  (A.8)  is  well 
known  in  the  literature  [2]. 

Next,  V^Vi^  may  be  identified  as  an  integrating  factor  in  equation  (3.7),  i.e. 

(A.9) 
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However,  from  (A.2),  (A.6)  we  obtain 


*”  %k ' 

Vj^j  =  {Aj  exp[iA:j(z  —  z)]  —  Bj  exp[— —  z)]) . 


Now,  substitution  of  (A.IO)  into  (A.9)  yields 

cjiW  =  * 


=  (|r)' 


v,Hz)ViHz) 

where  the  integrcil  is 

,  a2  d2  iexp[2i{kj  -  ki){z  -z)]-l)  ^  ,2  a  t>  (exp[2iA:j(2:  -  z)]  -  1) 

- ‘ - 2% - 

■  p2  .2  (exp[2i(fc;  -  kj){z  -  z)]  -  1)  .  „2  r.2  (exp[-2i(fcj  +  fcf)(z  -  z)]  -  1) 
^  '  2i{ki-kj)  -2i{kj+ki) 


(A.10) 

(A.11) 


242  (exp[2i(fcj  +  ki){z  -  z)]  -  1) 
2i{kj  +  ki) 

'j{z- 

2ik.i 


2iki 

z,A  T>  r>2  (exp[-2iA:i(z  -  ^)]  -  1)  , 

-2AjBjBi  - - +  AAjBjAiBi{z  -  z)  ]  . 


-2iki 


(A.12) 

When  j  =  I  m  (A.12),  indeterminate  quotients,  such  as  (exp[2i(A:j  —  ki){z  —  z)]  —  l)/2i{kj  —  ki), 
must  be  replaced  by  the  limiting  value  {z  —  z). 

Similarly,  is  an  integrating  frwjtor  in  (3.9).  We  obtain 

OiM  =  +  7"  l\vj(inv,‘Sfds'j  ,  (A.13) 

where  the  integral  can  be  evaluated  to  be 

(.4|(a;)2  ^)1  -  1) 

^2/D.^2  («^P[2i(A:j  +  kl){z  -  z)]  -  1)  ^  ^2  d*  (exp[2ifej(z  -  z)]  -  1) 


+AjiBtY 


2i{kj  +  At,*) 


2AjA*iBt 


2ikj 


,  r>2(  .*'>2  (exp[-2i(%  +  kf){z  -  Z)]  -  1)  2/  o*^2  “  kj){z  -  z)]  -  1) 

^  - 2i{kl  -  kj) - 

or2  a*p*  {exp[-2ikj{z  -  z)]  -  1)  „  .  p  /  4*^2  (exp[-2zA;f  (z  -  z)]  -  1) 

-2BjAiBi  - — - 2AjBj{A,)  - — - 


-2AjBjiBt) 


2  (exp[2ifc*(z  -z)]  -  1) 


2ikf 


+  4AiB. 


(A.14) 


Benjamin  White,  Exxon  Research  and  Engineering  Company,  Route  22E,  Annandale,  NJ  08801 
Werner  Kohler,  Department  of  Mathematics,  Virginia  Tech,  Blacksburg,  VA  24061 
Leonard  J.  Srnka,  Exxon  Exploration  Company,  P.  O.  Box  4778,  Houston,  TX  77210 


